The role of diffusion in branching and annihilation random walk models 
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Different branching and annihilating random walk models are investigated by cluster mean-field 
method and simulations in one and two dimensions. In case of the A -^ 2A, 2A — > model the 
cluster mean-field approximations show diffusion dependence in the phase diagram as was found 
recently by non-perturbative renormalization group method (L. Canet et al., cond-mat/0403423). 
The same type of survey for the A — > 2A, AA — > model results in a reentrant phase diagram, 
similar to that of 2A -^ 2,A, 4yl ^ model (G. Odor, PRE 69, 036112 (2004)). Simulations of 
the A — > 2A, 4:A — » model in one and two dimensions confirm the presence of both the directed 
percolation transitions at finite branching rates and the mean-field transition at zero branching rate. 
In two dimensions the directed percolation transition disappears for strong diffusion rates. These 
results disagree with the predictions of the perturbative renormalization group method. 



I. INTRODUCTION 

Phase transitions in nonequilibrium, dynamical sys- 
tems, which do not satisfy the detailed balance condition 
and do not possess hermitian Hamiltonian may appear in 
models of population, epidemics, catalysis, cooperative 
transport [1], enzyme biology [2] and markets for exam- 
ple [3] . Many of the known systems can be mapped onto 
some reaction-diffusion models, the behavior of them are 
the studied intensively in the past decades [1,4]. In these 
systems particle (A) creation, annihilation and diffusion 
processes compete and by tuning the control parameters 
phase transition may occur from an active steady state 
to an inactive, absorbing state of zero density. The sim- 
plest example of such models exhibiting phase transition 
are the branching and annihilating random walk models 
(BARW), in which offsprings are created by a single an- 
cestor: A ^ {n + 1)A and particles annihilate: 2A -^ 0. 

The classification of universality classes of nonequilib- 
rium systems is one of the most important tasks of statis- 
tical physics [5,6]. Universal scaling behavior may occur 
at continuous phase transitions like in equilibrium sys- 
tems and the corresponding n-point correlations are ho- 
mogeneous (but anisotropic) functions of space and time. 
In the past decades analytical and numerical studies ex- 
plored a large zoo of such classes [6] . One may pose the 
important question: which are the relevant factors de- 
termining these classes ? By inspecting the widespread 
literature it appears that besides the well known factors 
of homogeneous equilibrium models with short ranged 
interactions: i.e. spatial dimension, symmetries, conser- 
vation laws there are more things must be taken into ac- 
count. For example due to the possibility of transitions 
in low dimensions topological effects may play important 
role [7]. Furthermore initial condition - as a boundary 
in the time direction - can also be relevant [8,9]. 

The mean-field classes of general, 



site occupation number restricted models (with n > 1, 
m>l, k>0,l>Q and m~ I > 0) resulted in a series of 
different universality classes depending on n and m [10]. 
This means that above the upper critical dimension (dc) 
of the transition n and m are relevant parameters. In 
particular for the n — m symmetrical case the density 
of particles in the active phase, near the critical point 
(ctc > 0) scales as 
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with (3—1, while at the critical point it decays as 

p ex f-" , (3) 

with a = /3/j^ii — 1/n [11,10]. On the other hand for the 
n < m asymmetric case continuous phase transitions at 
zero branching rate CTc = occur characterized by 



(3 ~ l/{m — n), a = l/{m — 1) 



(4) 
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For n > m the mean-field solution provides first order 
transition. 

By going beyond the site mean-field approximation it 
turned out that the phase diagrams of the type eq.(l) 
models may contain other transition points with non- 
trivial scaling behavior. In previous papers [12,13] I in- 
vestigated the 2 A — > 3A, AA — > model by cluster mean- 
field approximations and simulations in one and two di- 
mensions and showed that the diffusion plays an impor- 
tant role: it introduces a different critical point besides 
the one obtained by site mean-field solution at CTc = 
(4). The non-trivial critical point at ac > 0, appearing 
by low diffusion rates exhibits the universal behavior of 
the transition of the 2A -> 3A, 2A ^ (PCPD) model 
owing to the generation of the effective 2A -^ reaction 
via the quick processes: 2 A -^ 3 A -^ AA -^ [14]. 

Very recent studies by Canet et al. [15,16] found simi- 
lar diffusion dependence in the A -^ 2A, 2A ^ BARW 
model. Using non-perturbative renormalization group 



method (NPRG) they obtained a different phase diagram 
than what was expected by perturbative renormahza- 
tion (RG) arguments [17]. Namely directed percolation 
(DP) transition were found in these models for d > 2 di- 
mensions with dc = 4. Furthermore those transitions at 
CTc > branching rate persist up to infinite dimensions 
provided X/D is greater than a threshold value. This 
threshold value was found to be zero for d < 3 and finite 
for d > 3. 

In this work I apply dynamical, cluster mean-field ap- 
proximations for that BARW model and by extrapola- 
tions to A^ ^ oo cluster sizes I determine the phase dia- 
gram in Id. Furthermore I explore the phase transitions 
of the A — > 2 A, AA -^ model with this method and by 
simulations in one and two dimensions. 
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FIG. 1. (a) Phase diagram of the A 
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determined hy N = 3,4, ..8 (right to left curves) cluster ap- 
proximations, (b) The {\c/D)tr critical endpoint values of 
the CTc > transitions as the function of 1/A'^. Boxes corre- 
spond to the A -^ 2A, 2A -^ 0, while circles to the A — > 2A, 
4:A -^ model. The line shows a fitting of the form (10) 



The cluster mean-field method applied for noncquilib- 
rium models first by [18,19] are based on the calculation 
of iV-block probabilities of the model. This method has 
successfully predicted the phase diagram of many sys- 
tems. I apply it for one dimensional, site restricted lat- 
tice versions of the models mentioned before. It is well 
known that such approximations predict the phase struc- 
ture qualitatively well in Id. It provides an efficient phase 
diagram exploration method and could give valid results 
in higher dimensions too (see example [13]). One can set 
up a master equations for the Pm block probabilities as 



dt 



= f {Pn{{s..})) 



(5) 



where the site variables take the values: Si = 0, A. Dur- 
ing the solution of these equations one estimates larger 
than N sized block probabilities by the maximum overlap 
approximation like 



^^+1(^1, ...Sjv+i) 



Pn{si, ...SAr)PAr(s2, ■■■SjV+l) 
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Taking into account spatial symmetries, for the maximal 
N = 8 approximations of this work one has to find the 
solution of equations of 136 independent variables. Using 
this method first I investigated the model defined by the 
transition probabilities 
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The steady state solutions have been determined for 
N — 1,2, 3..., 8 approximations and the corresponding 
steady state densities ps {a, D) are calculated numeri- 
cally. The phase transition points are obtained for sev- 
eral values and are plotted in the ad D vs. \c/ D pa- 
rameterization on plot (Fig. 1(a)) in order to be compa- 
rable with the results of [16]. The phase transition lines 
for A^ > 1 corroborate the phase diagram of [16], with 
{\c/D)th > threshold values and with linear shapes in 
the a/ D -^ 00 limit. By inspecting the numerical values 
of the {\c/D)th{N) thresholds values for TV = 2, 3. ..8 (see 
Table I) one can setup a hypothesis that in general 
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exactly. Therefore one can extrapolate that in the N -^ 
cxD limit {Xc/D)th{N) -^ (Fig. 1(b)) in agreement with 
the results of [16] for Id. So A^ > 1 cluster mean-field 
approximations, which can take into account the diffu- 
sion exhibit similar phase diagram as the one obtained 
by the NPRG method and in the asymptotic limit they 
even reproduce the zero threshold corresponding to Id. 
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FIG. 2. Steady state density results in the N — 7 approxi- 
mations of the A — > 2A, AA — > model. Different curves cor- 
respond to: D = l, 0.7, 0.6, 0.4, 0.371, 0.36, 0.2 (top-to bot- 
tom). The (Tc > critical point disappears for D* > 0.37083. 



The other system of interest in this study is the 



AA2A, 



4yl 
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model, for which perturbative RG predicts dc — 2/3 [17], 
therefore in all physical dimensions only type (4) mean- 
field transitions are expected at ac = 0. By solving the 



cluster mean-field equations up to iV = 8, with A = 1 — cr 
parameterization one obtains a reentrant phase structure 
for low diffusion rates as in case of the 2A — > 3 A, AA -^ 
model [12] (Fig. 2). I have determined the {Xc/D)th 
threshold values for each N and plotted them as the func- 
tion of 1/A'' (Fig. 1(b)). One can see a tendency towards 
level-off contrary to the results for the A -^ 2A, 2A -^ 
case. By extrapolating to the iV — > oo limit using a form 
similar, but more general than (8) 



{Xc/D)th{N) = a + b/iN-lY 



(10) 



the threshold does not vanishes : a = 1.2, b = 34.5, 
c = 2.14, therefore we may expect phase transitions with 
CTc > even in one dimension for any < D < 1 diffusion 
probabilities. 






FIG. 3. Simulation results for the steady state density 
of the A — > 2A, AA -^ model. Crosses correspond to 
D = 1.0, "-I-" signs to D = 0.2 diffusions in Id. Circles 
denote D — 0.01, squares D = 1 data in 2d. Error bars are 
smaller than the symbol sizes. The lines serve to guide the 
eye. The insert shows the data magnified in the neighbor- 
hood of the a — transition point. The solid line shows a 
power-law fitting with the exponent /3 = 0.33(1). 



To test the cluster mean-field results for the A -^ 2A, 
AA — > model I performed simulations in one dimension 
with lattice sizes: L — 10^ and with periodic boundary 
conditions. The simulations were started from homoge- 
neously filled lattices with probability 1/2. One elemen- 
tary Monte Carlo step consist of the following sub-steps. 
A particle and a direction are selected randomly. A parti- 
cle hopping is attempted in the given direction with prob- 
ability D provided the nearest neighbour (nn) site was 
empty. If the diffusion attempt was unsuccessful cither 4 
nn particles in the given direction are removed with prob- 
ability A or an offspring is created in an empty nn site in 
the direction selected before. The time (f) measured by 
Monte Carlo steps (MCS) is updated by At = 1/n, where 
n is the total number of particles. The density of particles 
p(t) is followed up to 10^ MCS (throughout the whole pa- 
per the time is measured by Monte Carlo steps (MCS)). 



The phase transition points have been located for several 
diffusion rates and are shown on Fig. 3. A qualitative 
agreement with the low-D results of the finite iV-cluster 
mean-field approximations can be observed. However the 
(Tc > transition persists even for D = 1 diffusion prob- 
ability in agreement with the N -^ oo extrapolations. 

The density decay near a CTc > critical point, at 
D = 0.2 has been investigated in more detail. The den- 
sity at (Tc is expected to decay as (3). Fig. 4(a) shows the 
local slopes of a defined as 



ae//(0 



lii[p{t)/pit/m)] 
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(where I used m = 2). As one can see curves with 
a > 0.8398 veer up (active phase), while those with 
a < 0.8398 veer down (absorbing phase). A clean scal- 
ing with a = 0.159(1) exponent can be observed at 
a = 0.83980(2) corroborating the 1-1-1 dimensional DP 
class exponent value (a = 0.159464(6) [20]). In the inac- 
tive phase p{t) vanishes exponentially in agreement with 
the DP behavior again. One may expect the same kind of 
transition all along the adD) > transition line. Indeed 
simulations showed that the density decays in a similar 
way at transition points with D = 0.01, 0.05, 0.09. 
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FIG. 4. (a) Local slopes of the density decay exponent a 
as the function of t^°'^ of the the one-dimensional A — > 2A, 
A A —> model at D = 0.2. Different curves correspond to 
a = 0.8399, 0.83985, 0.8398, 0.83975, 0.8397 (top to bottom). 
The critical point is located at ac = 0.83980(2). (b) The same 
as in (a) in two dimensions, for a = 0.1426, 0.14262, 0.14263, 
0.14264, 0.14265, 0.1427 (top to bottom). 

The simulations for this model were repeated in two 
dimensions in L = 4 x 10"^ linear sized systems with 
periodic boundary conditions. These simulations were 
started from fully occupied lattices. One elementary 
Monte Carlo step consists of the following processes. A 
particle and a number xi € (0, 1) are selected randomly; 
if xi < D a, site exchange is attempted with one of the 
randomly selected nn. The time is updated by 1/n. A 
particle and a number X2 G (0, 1) are selected randomly. 
li X2 < cr and if the number of nn empty sites is greater 
than 0, one new particle is created at an empty site se- 
lected randomly. If 2:2 > cr and the number of nn parti- 
cles is greater than 2, four randomly selected neighbor- 
ing particles are removed. The time (i) is updated by 



1/n again. The density of particles was followed up to 
tmax < lO"^ MCS. As one can see simulation data (Fig. 3) 
and the 7-point approximations (Fig. 2) fit qualitatively 
well (this is true for other A^ > 1 levels as well). In 
both cases for weak diffusion rates reentrant phase tran- 
sitions occur with Cc > 0, while for strong diffusions only 
a single phase transition at ctc = branching rate can be 
found. The density decay at D = 0.01 is analyzed near 
the phase transition point. The local slope figure shows 
a separatrix for the critical value CTc = 0.14263(1) as can 
be seen on Fig. 4(b). One can read-off the corresponding 
decay exponent a = 0.445(5), which agrees with the 2-t-l 
dimensional DP value {a = 0.4505(10) [21]). 

I also investigated the steady state behavior at the 
(Tc — transition. The steady state density in the active 
phase near the critical phase transition point is expected 
to scale as eq.(2). As the insert of Fig. 3 shows by apply- 
ing power-law fitting in the 10^^ < cr < 10^'^ region one 
obtains f3 ~ 0.33(1) both in one and two dimensions at 
different diffusion rates. This agrees with the mean-field 
value (4) for this model. The density at ctc = decays 
as p oc t~-^/^ trivially, dictated by the 4^ ^ process. 

In conclusion numerical evidence is provided that the 
A -^ 2A, AA -^ branching and annihilating random 
walk process exhibits diffusion dependent phase transi- 
tions contrary the expectations coming from perturbative 
rcnormalization. In particular A^-cluster mean-field ap- 
proximations (with extrapolations to the N ^ oo limit) 
resulted in a reentrant {D-a) phase diagram, with phase 
transitions at CTc > for < I? < 1 diffusion proba- 
bilities. Simulations have shown that along this tran- 
sition line DP critical behavior occurs in one and two 
dimensions. This type of critical behavior is the con- 
sequence of an effective A ^ reaction generated by 
A -^ 2 A -^ 3 A -^ 4 A ^ of slowly moving particles. 
In one dimension this line of phase transitions persist for 
all < D < 1 diffusions, while in 2d it disappears for 
high diffusion rates. For any D and d values a mean- field 
transition, characterized by /? = 1/3 occurs at CTc — 0. 

Similar reentrant phase diagram has been observed in 
case of the one dimensional A -^ 2 A, 3A -> model [22], 
in a variant of the NEKIM model [23] and in the 2A -^ 
3A, 4yl ^ model [12,13]. In all cases the diffusion 
competes with particle reaction processes, and the bare 
parameters should somehow form renormalized reaction 
rates which govern the evolution over long times and dis- 
tances. An other study using exact methods [24] showed 
that the particle density fluctuations undergo a diffusion 
dependent phase transition in the bosonic PCPD model 
for d > 2. A very recent, non-perturbative RG study 
[16] has found similar diffusion dependent phase diagram 
in the A -^ 2A, 2A -^ model. That work points out 
that non-perturbative effects arise and there is a thresh- 
old (X/D)th{d) above which DP, while below it a type 
(4) mean-field transition at CTc = appears. The present 
A^-cluster approximations confirm those results and even- 



tuate similar phase diagram as in [16] for any finite N, 
with a threshold value: {\c/D)th{N) = 2/{N - 1). In 
the N ^ oo limit this corroborates a vanishing threshold 
in Id [16]. 
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1 


0.666 


0.5 


0.4 


0.333 


0.286 



TABLE I. Numerical N-cIuster results for the threshold 
values of the A ^ 2A, 2A -* model 



